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Abstract 

Recently  there  has  been  much  interest  in  scaling  water  flow  and  species  transport  at  the 
continuum  level  to  the  watershed.  A  particularly  simple  and,  therefore,  appealing  approach  is 
based  on  the  water  mass  balance  at  the  hillslope  scale.  Such  models  require  parameterization 
of  closure  relations  (flux-storage  relations)  based  on  held  data.  In  several  recent  studies, 
this  data  was  instead  generated  by  steady-state  numerical  simulations  of  the  hillslope.  In 
this  work  we  focus  specifically  on  closure  relations  for  hillslope  water  balance  models  as  in 
previous  work,  but  we  use  transient  numerical  solutions  of  continuum-scale  models  of  the 
subdomain  to  generate  the  data.  Our  goal  is  to  study  the  effects  of  non-equilibrium  behavior 
on  time-averaged  flux-storage  relationships  for  the  hillslope.  We  show  that  for  simulated 
hillslopes,  the  flux-storage  relations  are  multivalued  for  a  wide  range  of  time  scales,  discuss 
how  this  situation  arises,  and  provide  some  alternative  parameterizations  of  the  flux-storage 
relations. 


1.  Introduction 

Continuum  mechanical  models  of  water  how  form  the  basis  of  the  standard  physical  de¬ 
scription  of  watershed  hydrology.  In  principle,  at  least,  the  physical,  mathematical,  and 
computational  components  of  a  continuum  scale  model  of  large  catchments  and  watershed 
have  been  available  for  three  decades  [20].  In  practice,  general  continuum  models  have  not 
been  capable  of  making  accurate  predictions  for  systems  at  the  scale  of  watersheds  [3].  Lim¬ 
itations  of  the  continuum  modeling  approach  include  theoretical  difficulties  in  consistently 
formulating  and  closing  continuum  models  for  all  significant  how  processes,  mathematical 
and  computational  difficulties  in  solving  the  resulting  three-dimensional  partial  differential 
equations,  and  perhaps  most  importantly  difficulty  obtaining  accurate  model  parameters  at 
an  appropriate  scale. 

Our  focus  in  this  work  is  on  the  simplified  representation  of  hillslope  subsurface  how  for  the 
purposes  of  watershed  modeling.  As  accurate  continuum  models  of  variably  saturated  how 
at  the  scale  of  hillslopes  are  computationally  demanding  and  require  extensive  characteri¬ 
zation  of  the  subsurface,  researchers  have  pursued  several  alternative  strategies  to  modeling 
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hillslopes.  One  approach  is  to  use  approximate  or  exact  analytical  solutions  of  continuum 
models  for  the  hillslope  based  on  simplifying  assumptions  such  as  soil  homogeneity  or  sim¬ 
plified  geometry  [39] .  Parameter  estimation  can  be  used  to  fit  such  models  to  real  hillslopes 
where  the  physical  parameters  in  the  model  are  then  interpreted  as  effective  parameters.  A 
yet  simpler  approach  is  to  parameterize  integral  mass  balance  equations  for  a  hillslope  region 
directly,  based  on  physical  data  or  numerical  simulations  [17,  18,  15,  31,  33]. 

We  will  review  some  of  the  predominant  watershed  modeling  approaches  in  more  detail 
in  section  2  to  provide  a  context  for  our  work  on  hillslope  models.  As  a  broad  range  of 
techniques  have  been  applied  to  the  watershed  modeling  problem,  we  direct  the  interested 
reader  to  several  recent  journal  issues  devoted  to  modeling  issues  at  the  watershed  scale 
[4,  41]. 

Briefly  stated,  our  objective  in  this  work  is  to  study  the  effect  of  system  transience  on 
hillslope  water  balance  models.  In  particular,  their  effect  on  the  flux-storage  relations  that 
such  models  require.  A  logical  approach  for  an  initial  study  is  to  use  detailed  numerical 
simulations  of  idealized  hillslopes,  in  which  case  the  data  for  the  water-balance  flux-storage 
relations  is  easily  obtainable  [15]. 

In  section  3  we  present  our  derivation  of  the  hillslope  water  balance  model,  which  combines 
ideas  from  [15]  and  [31].  In  section  4  we  present  our  model  of  the  hillslope  based  on  macro¬ 
scopic  continuum  equations  and  incorporating  the  geometry  of  the  hillslope  and  nonlinear 
submodels  of  unsaturated  flow  processes.  We  study  the  behavior  of  this  continuum  model  to 
guide  the  parameterization  of  the  exchange  terms  required  for  closure  of  the  hillslope  model. 
There  are  natural  limitations  to  using  continuum  models  to  obtain  closure  relations  for  the 
water  balance  models,  and  there  are  many  open  questions  about  how  to  account  properly 
for  capillary  forces  and  heterogeneous  soil  properties  at  an  appropriate  scale.  Nevertheless, 
the  continuum  approach  is  capable  of  approximating  the  dynamics  of  simple  systems  and 
should  provide  a  useful  benchmark  and  starting  point  for  watershed  scale  models.  Finally  in 
section  5  we  present  the  resulting  flux-storage  relationships  for  the  hillslope  water  balance 
model  based  on  fully  transient  macroscale  data. 

In  summary,  our  objectives  are 

(1)  to  construct  a  simulator  for  a  single  hillslope  based  on  macroscale  continuum  models 
of  subsurface  flow, 

(2)  to  obtain  flux-storage  relations  for  the  hillslope  using  data  from  the  continuum  model, 
and 

(3)  to  explore  the  feasibility  of  constructing  a  low-dimensional  model  from  these  relations. 

2.  Background 

The  long  term  objective  of  watershed-scale  hydrological  modeling  is  to  predict  the  response 
of  watersheds  to  input  of  precipitation  from  the  atmosphere,  extraction  of  moisture  via 
evapotranspiration  and  discharge  via  regional  groundwater  flow  and  channel/overland  flow. 
There  are  a  number  of  different  modeling  approaches.  First,  one  could  draw  on  the  large 
body  of  work  on  continuum  models  for  subsurface  and  surface  fluid  flow  and  transport 
models  to  obtain  a  coupled,  spatially  distributed  continuum  model  for  the  watershed  system. 
This  approach  was  outlined  and  implemented  as  far  back  as  1969  [20].  Second,  one  could 
formulate  a  finite-dimensional  model  of  the  watershed  using  various  techniques  to  yield  a 
model  appropriate  for  a  significantly  larger  scale  (possibly  in  time  and  space)  to  yield  a 
description  of  the  averaged  water  balance  [17].  We  have  in  mind  in  the  latter  case  simple 
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water  balance  models  as  well  as  complex  linkages  of  analytical  solutions  of  the  continuum 
equations  and  other  approximations  [2],  Over  the  last  decades  widely  used  models  have 
evolved  as  some  mixture  of  both  approaches,  but  we  nevertheless  divide  recent  and  historical 
approaches  into  two  categories:  continuum  models  and  process  models. 

The  continuum  modeling  approach  for  watersheds  is  based  on  coupling  three  dimensional 
continuum  model  equations  governing  all  relevant  processes  participating  in  watershed  hy¬ 
drodynamics.  Thus,  at  a  minimum  models  of  flow  in  porous  media  must  be  coupled  to  open 
channel  flow  and  boundary  conditions  reflecting  evapotranspiration  and  precipitation.  This 
approach  was  formulated  in  [20],  though  given  computational  and  mathematical  limitations 
at  the  time  the  subsurface  sub-models  were  one-  and  two-dimensional  approximations  to  the 
full  three-dimensional  system.  Advances  in  computing  power  and  numerical  analysis  have 
led  to  many  recent  attempts  at  modeling  catchments  and  hillslopes  using  three-dimensional 
subsurface  flow  models  most  notably  [29,  9,  27].  These  more  recent  models  address  some  of 
the  shortcomings  of  continuum  models  cited  by  [20]  by  incorporating  digital  elevation  data 
for  defining  the  spatial  domain,  and  employing  models  of  canopy  interception,  evapotranspi¬ 
ration,  and  soil  heterogeneity. 

In  spite  of  recent  advancements  in  continuum  modeling  of  the  various  watershed  flow  pro¬ 
cesses,  several  of  the  shortcomings  cited  in  [20]  and  in  more  recent  critiques  of  watershed 
modeling  [5,  7,  8,  3,  30]  remain.  We  break  down  these  shortcomings  into  two  groups:  1) 
the  inability  to  determine  the  physical  parameters  of  the  models  at  the  continuum  scale  for 
the  entire  watershed  from  either  in  situ  measurements  or  parameter  estimation,  assuming 
we  have  correct  physical  model  of  the  processes,  and  2)  inability  to  characterize  all  the  rel¬ 
evant  processes  with  a  rigorous  physical  model.  Both  barriers  to  progress  are  rooted  in  the 
variability  of  the  domain  and  boundary  conditions,  and  both  may  yield  partially  to  contin¬ 
ued  developments  in  modeling  and  measurement.  On  the  other  hand,  severe  limitations  on 
modeling  predictions  due  to  the  propagation  of  uncertainty  may  be  a  fundamental  limita¬ 
tion  for  continuum  modeling  at  this  scale.  Given  that  in  many  applications  finely  detailed 
knowledge  of  the  hydrodynamics  is  not  even  required  of  the  models,  many  researchers  have 
pursued  simplifications  of  the  full  three-dimensional  continuum  approach. 

Quasi-three-dimensional  models,  where  one  or  two  dimensions  have  been  integrated  to 
yield  a  simplified  continuum  model,  are  widely  used  for  describing  saturated  and  unsaturated 
flow  in  porous  media  [10,  11].  A  number  of  other  model  simplifications  still  strongly  tied 
to  the  continuum  approach,  including  an  array  of  upscaling  approaches  such  as  volume 
averaging  and  homogenization,  will  not  be  discussed  further  here.  The  simplifications  that 
we  have  termed  process-based  models  make  a  more  complete  break  with  continuum  models 
and  yield  directly  a  model  whose  solution  is  itself  finite  dimensional.  The  most  widely  used 
example  of  such  models  is  TOPMODEL  [2],  which  discretizes  the  watershed  based  on  a 
topographic  similarity  index  related  to  the  surface  slope  and  drainage  area.  TOPMODEL 
makes  use  of  analytical  solutions  of  continuum  equations  to  piece  together  a  description  for 
each  relevant  process. 

A  number  of  other  process-based  models  have  appeared  over  the  last  fifteen  years  [17,  18, 
38,  23,  44,  43,  12,  13,  24,  25,  26].  Agreement  between  process-based  models,  held  data,  and 
continuum  models  has  been  the  subject  of  a  several  studies  [40,  37,  36,  21],  The  results  of 
these  comparisons  have  been  mixed.  However,  given  that  the  much  simpler  process-based 
models  are  quite  good  at  predicting  at  least  a  subset  of  the  dynamics  recorded  in  the  held  and 
simulated  by  detailed  continuum  models,  improving  process-based  models  might  be  a  worthy 
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Figure  1.  A  watershed  partitioned  into  3  hillslopes.  Streams  are  shown  in 
bine  (solid  lines),  drainage  divides  in  red  (dotted  lines),  and  the  watershed 
boundary  in  black  (dashed  lines). 


avenue  of  research.  Significant  progress  has  been  made  on  analytical  approximations  to 
continuum-based  subsurface  flow  models  with  the  express  purpose  of  incorporating  dominant 
topographic  effects  [39]  into  process  models. 

Due  to  the  widespread  use  of  process-based  models  and  the  realistic  dynamics  produced  by 
continuum  models,  some  recent  research  has  focused  on  using  continuum  models  to  improve 
process  models.  Duffy  and  colleagues  studied  equilibrium  solutions  of  a  continuum  model  of 
unsaturated  flow  in  a  simplified  hillslope  [15,  16,  6]  in  order  to  extract  a  simplified  model  of  a 
single  hillslope  for  both  flow  and  species  transport.  Numerical  solutions  of  continuum  equa¬ 
tions  on  hillslopes  have  also  been  used  to  study  the  dependence  of  saturated  area  formation 
on  topographic  factors,  soil  properties,  and  rainfall  intensity  [28]. 

A  comprehensive  watershed  modeling  framework,  was  presented  in  [31,  33]  and  applied 
in  [32,  34],  This  formal  framework  for  finite  dimensional  models  of  watersheds  could  con¬ 
ceivably  unite  many  of  the  useful  features  of  process  models  into  a  more  rigorous  physical 
and  mathematical  theory  based  on  volume  and  time  averaging  as  well  as  thermodynamically 
constrained  closure  of  the  balance  equations.  Our  approach  will  take  the  time-averaging 
approach  used  in  the  framework  and  examples  presented  in  [31,  33,  32]  to  derive  models 
quite  similar  to  the  water  balance  models  in  [15,  16,  6]. 

3.  Hillslope  Water  Balance 

To  formulate  watershed  models  based  on  integrated  continuum  scale  quantities  we  must 
specify  an  integration  volume.  In  this  work  the  integration  volume  is  a  single  hillslope,  which 
is  a  subdomain  of  a  watershed  bounded  by  a  stream  reach,  the  land  surface,  and  a  set  of 
drainage  divides.  Any  watershed  can  be  decomposed  into  hillslopes  as  follows:  A  watershed 
has  an  associated  stream  network  consisting  of  all  streams  in  the  watershed  as  shown  for 
a  simple  watershed  composed  of  three  hillslopes  in  figure  1.  The  stream  network  can  be 
decomposed  into  segments  (channel  reaches).  Each  channel  reach  has  an  associated  drainage 
area  determined  by  the  topography,  and  each  drainage  area  can  itself  be  decomposed  into  a 
surface  (overland)  and  subsurface  flow  subdomains.  Thus  the  organization  inherent  in  the 
watershed  furnishes  a  decomposition  of  the  watershed  into  channel  reaches,  subsurface,  and 
overland  flow  [31]. 


5 


In  this  work  we  focus  simply  on  the  hillslope  water  balance,  in  particular,  on  the  flow  of 
water  out  of  the  hillslope  as  either  exfiltration  to  the  land  surface  or  base  flow  to  the  stream 
reach.  We  will,  however,  make  strong  limiting  assumptions  about  the  channel  reach  and 
land  surface  bounding  the  hillslope  in  order  to  constrain  the  scope  of  this  work  to  subsurface 
processes  as  was  done  in  [15].  These  assumptions  are  as  follows 

(1)  The  water  level  in  the  channel  is  constant. 

(2)  The  water  level  in  the  overland  flow  region  is  negligible. 

(3)  The  exchange  of  water  between  the  atmosphere  and  the  unsaturated  portion  of  the 
surface  is  determined  by  the  precipitation  rate. 

(4)  The  exchange  of  water  between  the  subsurface  subdomain  and  neighboring  subdo¬ 
mains  is  negligible. 

Thus  we  will  not  maintain  the  capability  to  simulate  any  feedback  to  the  subsurface  due  to 
significant  surface  ponding  or  water  level  fluctuations  in  the  channel,  nor  will  we  be  able  to 
simulate  so-called  infiltration  excess  overland  flow  or  regional  groundwater  flow.  While  the 
assumptions  are  likely  valid  for  some  hillslopes,  such  as  a  steep  hillslope  bordering  a  large 
channel  reach,  these  assumptions  are  only  invoked  to  allow  us  to  focus  on  the  subsurface 
dynamics  and  could  be  relaxed  if  the  continuum  scale  model  described  in  the  next  section 
was  coupled  to  continuum  models  for  open  channel  and  overland  flow.  In  the  context  of 
macroscale  continuum  models  of  the  subsurface  our  hillslope  is  bounded  by  a  low  permeability 
material  except  where  the  subdomain  intersects  the  channel  reach  and  the  land  surface. 


3.1.  Integration  in  Time.  Following  [31]  we  introduce  a  characteristic  time  scale  At.  The 
hillslope  water  balance  we  consider  will  be  a  time- integrated  mass  balance  calculated  in  terms 
of  continuum  variables  over  a  hillslope  with  volume,  V.s,  and  time  interval  T  =  [t  —  At,t  +  At\ 
of  length  T  =  2  At.  This  approach  is  different  from  that  used  in  [15]  where  only  steady-state 
solutions  of  the  underlying  continuum  model  were  considered,  and  allows  us  to  investigate 
the  effects  of  system  transience  on  our  ability  to  generate  unique  low  order  flux-storage 
relations. 


3.2.  Mass  Conservation  Equations.  Continuum  scale  quantities  for  water  in  the  subsur¬ 
face  are  density  p,  saturation  S,  and  porosity  e  where '  denotes  a  macroscale  quantity.  The 
time-averaged  water  mass  in  the  hillslope  is  then 

(1)  M——f  f  eSpdVdr 

T  Jt  Jvs 

where  Vs  is  the  subsurface  domain  of  the  hillslope.  We  now  factor  A4  into  physically  relevant 
quantities.  First  we  define  the  hillslope  volume,  porosity,  and  saturation 


(2) 

(3) 

(4) 


V  = 


e  = 


S  = 
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Note  that  with  this  definition  S'  =  1  if  and  only  if  S  =  1  throughout  the  entire  hillslope. 
The  hillslope  water  density  is 

(5)  p=W7sil^dvdT 

Thus  we  have  finally 

(6)  A4  =  pSeV  —  —  f  f  eSpdVdr 

T  Jt  JVs 

Note  that  if  e  is  constant  in  time  and  p  is  constant  in  time  and  space  then  e  and  p  =  p  are 
both  constants.  That  is,  if  the  medium  and  water  are  assumed  to  be  incompressible  then  the 
analogous  assumptions  hold  for  the  hillslope.  We  will  derive  mass  conservation  on  a  mass 
per  unit  volume  basis  so  we  define 

(7)  m  =  pSe 


If  there  are  no  internal  sources  of  mass  in  the  subdomain  then 

(8)  —  f  epSdV  +  f  n  •  epSvdA  =  0 

JVs  JAs 

where  As  is  the  boundary  of  Vs,  n  is  the  unit  outward  normal  on  As,  and  v  is  the  macroscopic 
velocity  of  the  water  phase  (i.e.  the  filtration  velocity  or  Darcy  velocity).  Taking  the  time 
average  and  dividing  by  the  subsurface  volume  V 

o)  ±  x  s  L  Mvdr +^vIrL&' MdAdT = 0 

We  can  interchange  differentiation  and  integration  to  obtain 

(1°)  MVdr  mdAdr  =  0 


or  simply 

(ii) 


dm 

~dt 


=  e 


A 


Using  the  variables  and  simplifying  assumptions  above  we  partition  the  exchange  term  eA 
into  terms  from  the  subsurface  region  (s)  into  the  overland  flow  (o),  channel  reach  (r),  and 
atmospheric  (a)  bounding  regions 


d  (pSe) 

(12)  7^  €so  T  6sr  T  &sa 

The  simple  hillslope  with  the  exchange  terms  labeled  is  shown  in  figure  2.  We  have  formulated 
a  hillslope  equation  for  mass  conservation,  which  now  needs  to  be  closed  by  finding  additional 
equations  describing  the  exchange  terms.  If  p  and  e  are  constant,  we  need  only  find  a  relation 
between  the  fluxes  e  and  the  storage  S.  Our  closure  approach  will  be  to  generate  data  using 
macroscale  numerical  simulations.  We  next  describe  the  macroscale  model. 
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Figure  2.  2D  Hillslope.  The  region  is  divided  into  stream  reach  (Vr)  and 
subsurface  (Vs)  subdomains.  The  water  table  divides  Vs  roughly  in  half,  and 
the  fluxes  out  of  each  boundary  segment  arc  labeled.  For  our  sign  convention, 
fluxes  out  of  the  domain  are  negative. 


4.  Continuum  Modeling  Approach 

We  will  follow  the  standard  macroscale  continuum  approach  for  modeling  flow  in  the 
subsurface  [1].  The  application  of  continuum  models  of  the  subsurface  to  hillslope  and 
watershed  systems  was  apparently  first  outlined  in  [20].  Notable  recent  efforts  are  described 
in  [42,  7,  29]. 

4.1.  Assumptions.  We  use  the  following  assumptions,  which  are  consistent  with  those  given 
above,  in  formulating  the  continuum  model  equations 

(1)  The  geometry  and  pressure  distribution  of  the  channel  is  constant,  in  particular  we 
assume  a  static  equilibrium  pressure  distribution. 

(2)  The  water  level  in  the  overland  flow  subdomain  is  negligible  and  constant,  in  partic¬ 
ular  the  pressure  is  atmospheric  for  the  thin  sheet  of  overland  flow. 

(3)  The  exchange  of  water  between  the  atmosphere  and  the  unsaturated  surface  of  the 
hillslope  is  determined  by  atmospheric  conditions. 

(4)  The  hillslope  is  symmetric  about  a  long  channel  reach.  We  will  consider  a  2D  slice 
of  hillslope  by  symmetry. 

(5)  Air  in  the  subsurface  is  at  constant  atmospheric  pressure  and  infinitely  mobile. 

4.2.  Boundary  and  Initial  Conditions.  The  assumptions  above  allow  us  to  ignore  simu¬ 
lating  overland  flow  and  channel  flow  and  to  simulate  the  hillslope  using  Richards’  equation 
[35].  The  effects  of  the  channel  and  overland  flow  subdomain  enter  only  through  the  bound¬ 
ary  conditions  that  we  use  for  Richards’  equation.  The  boundary  conditions  are 


(13) 

A _ 

s^-sr 

p  =  pg(ztop  -  z) 

(14) 

A _ 

•ST-SO 

v  •  n  =  max(k  ■  p  —  qp,  —qp) 

(15) 

■AsdV - 

v  •  n  =  0 

where  ztop  is  the  depth  of  the  channel  reach,  qp  is  the  precipitation  rate  and  A;  is  a  parameter 
representing  surface  conductivity  (lOm/d  in  this  work).  The  initial  conditions  were  static 
equilibrium  conditions  corresponding  to  qp  =  0.  We  use  a  mass  conservative  formulation 


Table  1.  Homogeneous  Hillslope  Parameters.  These  values  correspond  to  a 
sandy  sooil. 
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u 

3.010  x  10-1 

Sr 

2.799  x  10~2 

Ks  (m/day) 

5.040  x  10u 

a  (m-1) 

5.470  x  10u 

n 

4.264  x  10u 

Po  (kg/m3) 

9.982  x  102 

(3  (m  •  day2  /  kg) 

6.564  x  10“2U 

of  Richards’  equation  and  numerical  methods  that  are  described  in  [22]  along  with  closure 
relations  of  Mualem  and  Van  Genuchten.  The  mass  balance  and  closure  relations  are  given 


by 

(16) 

dpuS 

at  +V'M 

=  0 

(17) 

V 

=  —K  (Vp  —  pg) 

(18) 

P 

=  Poe* 

(19) 

S 

=  Se(l-Sr)  +  Sr 

(20) 

Se 

=  [1  +  (a  max(—  p,  0)n]-m 

(21) 

K 

=  Ksy/se  [i  -  (i  -  symr] 

where  po,  (3,  u,  g,  Ks,  Sr,  a,  n  and  m  =  1  —  1/n  are  constants. 


4.3.  Hillslope  Simulator.  For  this  work  the  hillslope  domain  Vs  is  a  10m  x  100m  rect¬ 
angular  region  with  the  long  side  tilted  at  an  angle  7t/4  with  the  horizontal.  The  channel 
reach  boundary  is  applied  at  the  left  hand  side  so  that  ztop  =  10  sin  7t/4.  The  fluid  and  soil 
parameters,  which  correspond  to  a  homogeneous  sand,  are  given  in  table  1.  For  complete¬ 
ness  we  also  tested  slopes  of  7t/5  and  7t/6,  correlated  random  fields  for  the  soil  parameters 
based  on  the  Miller-similar  approach  with  a  variance  of  1.0,  and  a  block  heterogeneous  slope 
described  in  detail  in  [19] .  Since  the  features  we  describe  in  what  follows  are  apparent  in  the 
simple  homogeneous  slope  described  above  and  across  the  range  of  hillslopes  we  studied,  we 
will  henceforth  deal  only  with  the  simple  homogeneous  slope  with  slope  7t/4. 

5.  Watershed- Scale  Closure  Relations 

The  simplest  approach  to  closing  the  equation  for  the  subsurface  would  be  to  assume  a 
constant  watershed-scale  density  and  porosity  and  then  determine  the  exchange  terms  as 
functions  of  the  watershed-scale  saturation.  For  instance,  given  suitably  defined  watershed 
scale  pressure  p  we  might  assume  that  a  “watershed-scale  Darcy’s  law”  holds  [33]  and  that 
furthermore  p  can  be  determined  from  S.  Since  the  pressure  in  the  channel  reach  is  constant 
this  would  yield  the  parameterization 

(22)  esr  =  esr(S) 

In  order  to  test  the  hypothesis  that  esr  =  esr(S)  for  our  simple  system  we  ran  a  series 
of  simulations  to  equilibrium  starting  from  both  fully  saturated  and  fully  dry  (equilibrium 
for  P=0.0)  initial  conditions  for  precipitation  rates  P  =  0.0,  0.05,  0.1, ...,  1.0  (m/day).  We 


Figure  3.  esr  vs.  S  for  a  watershed  time  scale  of  T  —  3.2  days,  6  =  7t/4, 
homogeneous  slope.  Equilibrium  values  are  denoted  by  *.  Each  line  represents 
drainage  and  infiltration  to  equilibrium  under  constant  precipitation. 


collected  the  spatially  integrated  fluxes  esr  and  eso  as  well  as  watershed  scale  saturation  S  at 
0.1  day  increments  for  100  days,  which  was  approximately  enough  time  for  the  hillslope  to 
reach  equilibrium  with  the  precipitation  from  both  sets  of  initial  conditions.  To  obtain  the 
watershed-scale  variables  we  approximated  time  integrals  of  these  quantities  for  several  time 
scales,  T  =  0.1,  0.2, 0.4, . . . ,  102.4  days,  using  the  midpoint  rule.  Note  that  equilibrium  values 
are  not  affected  by  the  temporal  upscaling  since  the  system  is  assumed  to  be  at  equilibrium  for 
t  >  100  days.  That  is,  the  equilibrium  state  dominates  the  time-averaging  for  T  large  enough. 
In  order  to  preserve  this  property  in  our  discrete  integral  approximations  it  is  necessary  to 
add  sufficiently  many  copies  of  the  equilibrium  value  (T  days  worth  of  equilibrium  values) 
of  each  variable  to  the  100  day  sequence.  A  plot  of  esr  versus  S  is  given  in  figure  3  for 
a  watershed  time  scale  of  T  =  3.2  days.  The  result  is  that  esr  cannot  be  parameterized 
simply  as  a  function  of  the  watershed  saturation  since  esr(S)  is  multivalued.  If  we  were  to 
use  only  equilibrium  data  (an  asterisk  denotes  the  equilibrium  value  of  esr  in  figure  3),  our 
data  reproduces  the  result  in  [15]  where  the  exchange  term  was  parameterized  as  a  low  order 
polynomial  of  the  saturated  storage.  Furthermore,  it  was  noted  in  [15]  that  for  moderate 
storage  the  equilibrium  flux  to  the  channel  reach  is  roughly  linear  in  the  saturated  storage, 
and  our  results  demonstrate  the  same  behavior  if  we  partition  the  hillslope  storage  S  into 
unsaturated,  Su,  and  saturated,  Ss,  storage  as  in  [15].  In  our  model  a  low  order  polynomial 
in  the  hillslope  saturation  S  would  likely  suffice  for  a  moderate  range  of  precipitation  rates. 
Regardless  of  whether  the  hillslope  is  modeled  as  a  single-state  or  two-state  model,  however, 
the  flux  terms  remain  multivalued.  In  the  latter  case  the  ( Su ,  Ss,  esr)  data  form  a  multivalued 
surface  (a  surface  that  overturns).  As  should  be  clear  from  the  fact  that  the  equilibrium 
states  dominate  the  time  integration  for  large  enough  time  scales,  increasing  the  watershed 
time-scale  T  tends  to  collapse  the  multivalued  flux  functions  onto  the  equilibrium  curve. 
The  source  of  the  non-uniqueness  above  in  fluxes  with  respect  to  hillslope  storage  is  that 
disturbances  in  pressure  propagate  with  infinite  speed  due  to  the  non-degenerate  parabolic 
form  of  the  governing  equations  (Richards’  equation),  and, therefore,  for  any  given  watershed 


Figure  4.  esr  vs.  (S,  P )  for  T  =  3.2  days.  The  height  of  the  surface  gives  flow 
into  the  subsurface  from  the  stream.  The  equilibrium  contour  is  superimposed 
in  blue.  The  surface  is  single  valued  for  all  ( S,P ). 
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Figure  5.  eso  +  esa  vs.  (. S,P )  for  T  =  3.2  days.  The  height  of  the  surface 
gives  flow  into  the  subsurface  from  the  surface,  including  both  atmospheric 
fluxes  into  the  domain  and  exfiltration  from  the  saturated  subsurface.  The 
equilibrium  contour  is  superimposed  in  blue. 


scale  saturation  S  (or  likewise  given  values  of  Ss  and  Su)  infinitely  many  values  of  the  fluxes 
at  the  r  and  o  boundaries  can  be  generated  simply  by  varying  the  atmospheric  flux  esa. 

One  route  to  parameterizing  esr  is  then  to  add  another  dependency  that  reflects  the 
boundary  conditions.  Figure  4  plots  the  same  data  versus  both  S  and  the  precipitation  rate 
P,  which  shows  a  single  valued  surface  for  esr.  Figure  5  gives  the  total  flux  at  the  hillslope 
surface  as  a  function  of  ( S ,  P)  as  well. 
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Figure  6.  esr  vs.  (S,  P )  for  T  =  3.2  days  .  The  baseflow  driven  by  the 
dynamic  precipitation  time  series  is  shown  in  purple. The  fluxes  often  do  not 
lie  on  the  surface,  particularly  during  drainage. 

3.2  Day  Average 


5.1.  Realistic  Atmospheric  Forcing.  The  data  set  for  figures  3,  4,  and  5  consists  of 
monotonic  drainage  and  wetting  experiments,  and  a  more  realistic  case  would  be  relatively 
short  precipitation  events  interspersed  with  long  drying  periods.  To  this  end  we  obtained 
a  year  of  precipitation  data  from  a  site  in  North  Carolina  consisting  of  15  minute  cu¬ 
mulative  precipitation  measurements  (Station:  Clinton  2  NE,  COOP:  311881,  year  1999, 
http://www.ncdc.noaa.gov).  If  we  use  this  data  to  drive  the  hillslope  then  the  flux  becomes 
significantly  more  complicated.  We  plot  the  values  from  the  variable  precipitation  run  as  the 
purple  trajectory  superimposed  on  the  the  previous  data  in  figure  6  and  7.  The  exchange 
term  is  yet  again  multivalued  for  both  the  simple  storage  parameterization  as  well  as  the 
more  complex  parameterization  of  flux  with  respect  to  (. S,P ).  Either  approach  would  be 
biased  under  some  conditions  toward  more  base  flow  than  actually  occurs.  The  amount 
of  over-prediction  relaxes  to  zero  over  time.  Choosing  larger  watershed  time  scales  damps 
some  of  this  behavior,  but  even  for  the  largest  time  scale  in  this  study  (T  =  102.5)  days  the 
flux-storage  behavior  is  still  quite  complex  (figure  8).  More  importantly,  the  time-integrated 
flux-storage  data  for  the  hillslope  under  dynamic  forcing  does  not  collapse  onto  the  equilib¬ 
rium  curve.  In  other  words,  the  average  state  of  the  dynamically  forced  hillslope  is  not  an 
equilibrium  state,  which  is  well  known  (c.f.  [14]). 

5.2.  Watershed  Model  Closure  and  Simulations.  To  eliminate  as  many  sources  of  error 
as  possible  we  simply  used  piecewise  bilinear  splines  to  generate  closure  relations  from  the 
flux-(S',  P )  surface  data.  For  simplicity  we  lumped  the  surface  fluxes  eso  and  esa  into  a 
single  term,  eso  =  esr  +  esa  (figure  5).  Lastly,  to  compensate  for  the  fact,  discussed  in  the 
previous  section,  that  the  dynamically  forced  hillslope  drains  more  slowly  than  the  numerical 
experiments  used  to  generate  the  flux-storage  data,  we  shifted  the  input  to  the  bilinear  spline 
for  esr  by  a  constant,  S*.  In  summary,  the  watershed-scale  model  for  the  simple  hillslope  is 

dS 
dt 


(23) 


e„(S  -  S')  +  e,a(S) 
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Figure  7.  esr  vs.  (S)  for  T  =  3.2  days;  The  surface  driven  by  the  dynamic 
precipitation  time  series  is  shown  in  purple.  Over  time,  the  base  flow  appears 
to  relax  to  the  trajectory  corresponding  to  zero  precipitation. 


Figure  8.  esr  vs.  (S)  for  T  =  102.5  days;  The  response  driven  by  the  dynamic 
precipitation  time  series  is  shown  in  blue.  For  large  time  scales  the  dynamic 
data  moves  away  from  the  zero  precipitation  trajectory. 


For  this  work  we  simply  set  S*  =  me  an  (5'  —  S)  where  S  is  the  output  of  the  model  with 
S*  =  0,  and  S  is  the  data  generated  by  the  continuum  model.  We  solved  the  ordinary 
differential  equation  above  using  the  MATLAB  routine  ode23  with  a  maximum  time  step  of 
min(T/2, 1)  days.  The  output  of  the  model  and  the  data  generated  by  the  continuum  model 
are  presented  for  two  time-scales,  T  =  0.1  and  T  =  6.4  days,  in  figures  9  and  10.  In  table  2 
we  present  a  measure  of  error,  ||S'i(t)  —  -S' (t)  || oo ,  where  is  the  output  of  the  watershed-scale 
model  as  well  as  the  parameter  S*. 


Figure  9.  Model  comparison,  T  —  0.1  days.  The  over-prediction  of  base  flow 
is  particularly  apparent  after  large  precipitation  events. 
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Figure  10.  Model  comparison,  T  =  6.4  days.  The  time  averaging  damps  out 
small  errors  but  has  no  significant  effect  on  errors  due  to  large  precipitation 
events  and  long  dry  periods. 


Table  2.  Watershed-scale  model  errors  and  coefficient, S'* 


T 

0.1 

0.2 

0.4 

0.8 

1.6 

3.2 

0.0474 

0.0469 

0.0461 

0.0448 

0.0409 

0.0325 

s* 

0.0146 

0.0146 

0.0146 

0.0145 

0.0143 

0.0138 

T 

6.4 

12.8 

25.6 

51.2 

102.4 

||5f  —  5|  lnf 

0.0174 

0.0173 

0.0240 

0.0231 

0.0285 

s* 

0.0137 

0.0150 

0.0246 

0.0333 

0.0410 
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6.  Discussion 

In  response  to  a  change  in  forcing  conditions,  two  processes  occur  in  the  subsurface  1)  the 
redistribution  of  moisture  and  2)  the  propagation  of  pressure  disturbances.  The  former  is 
associated  with  a  very  slow  time  scale-hours  to  months  (depending  on  soil  properties,  initial 
and  boundary  conditions),  while  the  latter  occurs  nearly  instantaneously.  As  modeled  by 
Richards’  equation,  the  hillslope  boundary  fluxes  are  strongly  dependent  on  the  pressure 
gradient  and  the  moisture  distribution  clue  to  the  highly  nonlinear  form  of  the  macroscale 
closure  relations.  The  fast  time  scale  of  pressure  signals  produces  a  kind  of  non-locality  or 
rate-dependence  in  the  behavior  of  the  watershed-scale  system;  the  behavior  of  the  subsurface 
as  a  buffer  in  the  hydrology  of  the  watershed  depends  not  only  on  the  water  stored  in  the 
subsurface  but  also  on  the  rate  at  which  water  is  being  supplied  to  the  subsurface  since 
that  rate  affects  the  pressure  gradient  across  the  hillslope.  The  rate-dependent  effect  can  be 
incorporated  into  the  water  balance  model  by  parameterizing  the  fluxes  as  functions  of  both 
S  and  P.  The  slow  (and  variable)  time  scale  of  moisture  redistribution  produces  a  lag  or 
memory  effect  that  cannot  be  fully  quantified  with  the  simple  approaches  we  investigated. 
This  effect  is  particularly  apparent  under  drying  conditions. 

7.  Conclusions 

We  investigated  several  approaches  to  parameterizing  flux-storage  closure  relations  for  a 
time-averaged,  integrated  hillslope  water  balance.  The  approaches  are  not  able  to  reproduce 
all  of  the  complex  dynamics  of  a  simulated  hillslope  driven  by  a  year  long  sequence  of 
natural  precipitation.  The  complex  dynamics  for  the  simulated  hillslope  derive  from  the 
nonlinear  parabolic  form  of  the  governing  equations  at  the  continuum  scale.  While  some 
dynamic  effects  related  to  variable  precipitation  are  incorporated,  those  particularly  related 
to  the  slow  redistribution  of  moisture  within  the  hillslope  are  not.  The  time-averaging  over 
large  time  scales  may  be  useful  in  recovering  some  accuracy  if  only  long  term  averages  are 
required  of  integrated  hillslope  model.  The  hillslope  simulator  in  this  work  was  quite  simple. 
Phenomena  in  natural  hillslopes  such  as  hysteresis,  heterogeneity,  macropores,  and  fractures, 
could  conceivably  dampen  or  exacerbate  some  of  the  effects  we  noted.  Further  work  is  needed 
in  constructing  more  complex  continuum  models  of  the  hillslopes  as  well  as  in  formulating 
more  rigorous  upscaling  techniques  for  this  problem. 
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